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PREFACE 


The  work  described  here  was  performed  under  the  program  for 
analytical  crashworthiness  prediction  at  the  Transportation 
Systems  Center  sponsored  by  the  Research  Institute  of  the 
National  Highway  Traffic  Administration  of  the  U.S.  Depart- 
ment of  Transportation.  This  program  is  intended  to  provide 
engineering  data  that  can  be  applied  to  establishing  regula- 
tions relating  to  vehicle  collision  performance  to  improve 
motor  vehicle  safety.  The  principal  author,  Dr.  J.  Rossettos, 
Associate  Professor  of  Mechanical  Engineering  at  Northeastern 
University,  held  a temporary  appointment  as  a staff  consultant 
to  the  Transportation  Systems  Center  during  the  summer  of  1973. 
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1.  INTRODUCTION 


Structural  crashworthiness  plays  an  important  role  in  the  NHTSA 
mission.  There  are  definite  requirements  for  analytic  models  which 
can  serve  as  interpolation  tools  in  conjunction  with  crash  testing. 
These  simulation  models  can  also  serve  as  design  feasibility  and 
evaluation  tools. 

The  Transportation  Systems  Center  has  been  asked  by  NHTSA 
to  provide  support  in  the  formulation  and  implementation  of  such 
analytic  models  and  computer  programs  in  order  to  predict  vehicle 
crashworthiness.  This  support  will  include  the  acquisition,  im- 
plementation and  extension  of  existing  computer  codes. 

Some  promising  analytic  tools  are  the  Calspan- Shieh  two-dimen- 
sional frame  analysis  program,  the  Battel le-FMCCM  lumped  mass  pro- 
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gram,  and  the  Lockeed  three-dimensional  KRASH  program. 

In  recent  years  various  simulation  programs  have  been  developed 
to  model  the  dynamic  structural  response  under  vehicle  impact  con- 
ditions. The  models  vary  from  the  very  simple  which  can  give  only 
average  features  of  the  overall  response,  to  the  rather  complex, 
where  greater  detail  in  the  response  can  be  provided.  The  simplest 
form  of  analytic  simulation  to  date  has  been  embodied  in  simplified 
spring-mass  models  with  2-3  lumped  masses  and  less  than  ten 
degrees  of  freedom,  while  generalized  resistances  are  made  to 
represent  gross  vehicle  structural  properties.  An  example  of 
such  models  is  given  in  Reference  4.  Correlation  with  tests  de- 
pends heavily  on  making  a judicious  choice  for  the  parameters  which 
measure  the  generalized  resistances.  Therefore,  their  use  as  pre- 
dictive tools  is  limited  although  they  can  be  used  to  establish 
general  behavior.  A good  example  of  the  proper  use  of  such  models 
is  given  in  Reference  5. 

The  next  step  beyond  the  simplified  spring-mass  model  exists 

2 

in  the  BCL  simulation  program.  Four  masses  are  represented 
together  with  35  individual  nonlinear  resistances.  There  is  a 
restriction  to  unidirectional  motion.  Judgment  in  the  selection 
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of  mass  and  resistance  parameters  is  required,  especially  in 
interpreting  the  crush  data  to  be  used  for  the  six  different 
types  of  force-deformation  curves  which  are  available.  Of 
course,  this  offers  a degree  of  flexibility  to  the  user. 

There  are  certain  simulation  programs  which  are  labelled 
as  "hybrid"  models  because  they  require  as  necessary  input,  ex- 
perimental crush  data. ^ Correlation  of  the  static  deformation 
mode  with  the  dynamic  mode  is  at  present  a very  difficult  problem, 
so  that  extrapolation  to  other  environments  is  not  assured,  and 
again  experience  and  judgment  are  important  for  any  reasonable 
prediction  capability. 

12  3 7 

The  next  step  in  complexity  involves  the  frame  models  ’ ’ ’ 
which  are  comprised  of  a large  number  of  beam  elements  and  lumped 
masses.  Increasing  the  number  of  degrees  of  freedom  would,  of 
course  provide  increased  ability  for  the  evaluation  of  detailed 
structural  response  of  components  in  vehicle  impacts.  The  simplest 
of  the  frame  models  is  the  Calspan-Shieh  program,'*'  a two-dimensional 
model  which  provides  for  elasto-plastic  response  by  the  use  of 
plastic  hinges.  The  plastic  hinge  idea  is,  of  course,  a simplified 
approach  to  yielding  of  a beam  cross  section  since  details  of  the 
stress  distribution  over  the  cross  section  are  not  taken  into  ac- 
count. Correlation  with  limited  tests  has  been  shown  to  be  satis- 

3 

factory.  Another  frame  model,  KRASH  program,  was  originally 
developed  for  aircraft  structure,  and  in  principle,  it  can  be  made 
to  apply  to  vehicle  impact.  It  can  be  regarded  as  an  extension 
of  the  BCL  model,  consisting  of  lumped  masses  connected  by 
straight  beam  elements,  where  each  mass  has  three  translational 
and  three  rotational  degrees  of  freedom.  The  codes  can  include 
energy  absorber  devices  and  seat  collapse  mechanisms.  The  large 
deformation  characteristics  are  treated  by  piecewise  linearization, 
whereby  the  linear  stiffness  matrix  is  adjusted  for  plasticity  at 
each  time  step  by  multiplying  by  a stiffness  reduction  factor. 

This  factor  is  determined  from  static  crush  data,  so  that  again 
experimental  difficulties  similar  to  the  Kamal  model^  exist,  since 
static  and  dynamic  mode  behavior  is  not  easily  correlated. 
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A more  general  element  frame  model  is  the  three-dimensional 
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CRASH  program,  which  contains  additional  features  not  found  in 
the  last  two  models.  No  prior  assumption  on  the  locations  of 
plastic  hinges  is  required.  Moments  and  forces  at  the  nodes 
are  computed  by  numerical  integration  of  the  stress  distribution 
over  the  cross  section.  The  actual  stress  - strain  behavior  of  the 
material  may  be  used  directly  but  the  stress  state  must  be  monitored 
at  various  locations  accross  the  cross  section.  It  is  clear  that 
the  additional  effort  and  computer  time  required  for  this  program, 
still  does  not  allow  prediction  of  detailed  vehicle  response  much 
beyond  the  capability  of  the  previous  two  models.  This  is  because 
the  frame  concept  cannot  really  model  an  entire  vehicle  body,  and 
also  since  local  deformation  of  the  cross  section  and  joint  inef- 
ficiency are  not  accounted  for. 

TSC  and  the  University  of  Michigan  under  NHTSA  contract  con- 
clude that  there  is  a need  for  a hybrid  finite  element  program, 
which  would,  for  instance,  incorporate  shell,  frame,  lumped 
parameter  and  finite  difference  models.  These  models  would  in 
fact  form  modules  for  the  overall  program,  where  each  module  can 
be  regarded  as  a "super  - element" . It  should  be  pointed  out  that 
the  use  of  various  simulation  programs  will  be  dictated  by  the 
particular  impact  situation  to  be  modelled.  For  instance,  in  some 
low  speed  impact  situations  where  the  bumpers  alone  may  be  involved, 
the  simple  BCL  model,  or  in  combination  with  the  Cal  span - Shieh 
model  may  be  sufficient  to  handle  the  significant  features.  In 
any  case,  it  is  clear  that  the  three-dimensional  frame  structure 
will  form  an  important  module,  and  it  will  be  the  main  concern  of 
the  present  report.  Figure  1-1  shows  a typical  vehicle  where  some 
of  the  module  representations  are  indicated.  A spring  suspended 
large  mass  which  may  be  the  engine  is  shown,  together  with  sup- 
porting frame  and  plate  and  shell  portions. 

TSC  has  become  familiar  with  the  Calspan  model  and  has 
simplified  the  input  to  the  two-dimensional  program,  so  that 
it  can  be  run  by  personnel  with  little  knowledge  of  the  de- 
velopment details.  The  present  report  summarizes  the  analytic 
requirements  for  extension  of  the  Calspan  two-dimensional  frame 
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to  a three-dimensional  frame.  A future  document  will  discuss  the 
programming  requirements  for  implementing  the  analytical  develop- 
ment. Program  modifications  made  by  TSC  to  expedite  usage  by 
engineers  will  also  be  described. 


Figure  1-1.  Hybrid  Vehicle  Model 
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2,  ANALYSIS  OF  THE  THREE-DIMENSIONAL  FRAME 


The  three-dimensional  frame  is  a necessary  module,  since  it 
is  possible  with  it  to  simulate  at  least  the  significant  features 
of  the  large  dynamic  plastic  deformation  and  geometry  changes  of 
the  vehicle.  This  is  done  by  means  of  an  appropriate  breakdown 
into  straight  beam  elements  which  are  connected  at  nodes,  with  in- 
ertial effects  treated  by  means  of  lumped  point  or  rigid  body 
masses  at  the  nodes. 

With  a sufficient  number  of  beam  elements  and  nodes,  one 
should  be  able  to  obtain  in  adequate  detail,  displacement  and 
acceleration  time  histories  of  nodal  positions  which  represent 
relatively  important  points  in  the  vehicle  structure,  including 
of  course  the  three-dimensional  motions  of  the  occupant  compart- 
ments. By  combining  this  information  with  up-to-date  biomechanics 
data,  one  can  evolve  an  estimate  of  crashworthiness  in  a particular 
environment . 

By  use  of  the  frame  model,  the  basic  overall  structural 
dynamic  response  phenomenon  can  be  identified.  Therefore  with 
some  experience  and  careful  interpretation,  it  can  also  be  used 
in  a very  important  way  to  establish  the  usefulness  of  any 
selected  experimental  parameters  in  taking  and  using  test  data, 
so  that  costs  of  tests  are  minimized. 

In  order  to  define  the  various  quantities  to  be  used  in  the 
analysis  a hypothetical  frame  is  shown  in  Figure  2.1.  It  should 
be  clear  that  in  this  figure  not  all  nodes  and  elements  are  neces- 
sarily shown.  This  is  to  avoid  cluttering  the  figure,  since  its 
purpose  is  mainly  to  define  nomenclature. 

With  reference  to  Figure  2-1  the  following  definitions  are 

used : 

Nodes  - are  fictitious  bodies  to  which  two  or  more  beam 
elements  may  be  connected.  The  end  of  each  element  con- 
nected to  a particular  node  takes  on  the  displacements  and 
rotations  of  that  particular  node.  Plastic  hinges  and 


5 


Figure  2-1.  Frame  Geometry  and  Notation 
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external  forces  may  be  concentrated  at  nodes.  Lumped  masses 
may  or  may  not  be  assigned  to  a node.  For  instance,  the 
locations  numbered  1,  2,  3,  ....  24,  are  nodes.  The  solid 
circles  have  assigned  lumped  masses  while  the  hollow  circles 
do  not.  In  a three-dimensional  analysis  there  are  six  degrees 
of  freedom  per  node,  so  that  for  n nodes  a system  of  6 n 
equations  would  need  to  be  solved. 


Beam  elements  - are  uniform  straight  members  connecting  two 
nodes.  In  Figure  2-1  beam  elements  are  denoted  by  the 
circled  numbers  (i.e.,  ©.©.® )• 


Beam  members  - are  actual  beams  which  make 
They  differ  from  beam  elements  in  that  two 
can  occur  on  a beam  member.  For  instance, 
the  lines  between  nodes  2 and  12  and  nodes 
beam  members. 


up  the  frame, 
or  more  nodes 
in  Figure  2-1, 
12  and  16  are 


The  physical  assumptions  of  our  model  are  essentially  those  adopted 
in  the  CALSPAN  model  of  Reference  1.  The  implications  of  these 
assumptions  are  as  follows.  All  plastic  deformation  is  to  occur 
at  the  nodes  of  our  beam  element.  The  location  of  hinges  must  be 
chosen  a priori,  and  this  gives  the  lengths  of  our  elements.  The 
plastic  hinge  is  operative  when  appropriate  stress  resultants  and 
bending  moments  lie  on  a given  yield  surface  for  the  cross  section. 
Perfectly  plastic  behavior  is  assumed  so  that  material  strain 
hardening  is  neglected.  Also,  the  stress  resultants  at  the  cross 
section  at  initial  yield  are  not  significantly  different  from  those 
at  the  fully  plastic  section.  Finally,  the  frame  structure  may 
undergo  large  rigid  body  translations  and  rotations.  These 
assumptions  are  felt  to  be  reasonable  for  mild  steel,  thin 
members,  and  vehicle  type  frame  loadings. 

In  the  analysis,  the  matrix  displacement  method  of  frame 
analysis  is  to  be  used,  and  the  dynamic  problem  is  reduced  to  a 
nonlinear  initial  value  problem,  which  is  governed  by  simultaneous 
second  order  differential  equations  in  time.  The  solution  is  car- 
ried out  incrementally,  and  the  coefficients  of  the  6 n dependent 
variables  (n  = number  of  nodes  in  the  frame)  are  updated  at  each 
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time  increment  to  take  into  account  the  possible  initiation  of 
either  a new  plastic  loading  or  unloading  condition  at  some  of 
the  nodes.  This  means  that  the  stiffness  matrix,  which  will  be 
defined  shortly,  is  to  be  modified  at  each  time  step. 

In  the  next  few  paragraphs,  we  will  define  the  coefficients 
which  enter  into  the  equations  of  motion  and  the  dependent 
variables  whos.e  solution  is  sought.  Since  the  basic  building 
block  of  our  model  is  the  beam  element  let  us  first  describe 
the  forces  acting  on  it  at  the  nodes,  and  the  resulting  deforma- 
tions and  displacements  which  it  experiences. 

If  we  imbed  a local  coordinate  system  at  each  end  face  (or  • 
node)  of  the  element,  we  can  describe  large  rotations  (rigid  body 
and  deformation  type)  of  the  beam  elements  and  their  nodes  by 
studying  how  such  local  coordinate  systems  translate  and  rotate 
with  respect  to  a fixed  (global)  system.  In  regard  to  Figure  2-2 
line  1 joins  the  beam  element  end  points  (or  nodes)  at  faces  "a" 
and  "b".  The  global  coordinates  of  nodes  a and  b are  X^  , , Z^ 

and  X2 , Y?,  Z2,  respectively.  The  quantities  N and  M are  stress 
resultant  and  moment  vectors.  Note  that  from  here  on  a wavy 
symbol  under  a letter  will  denote  a vector  or  a matrix,  and  will 
be  clear  from  the  context.  In  Figure  2-2  (b)  1 is  a unit  vector 

~3 

along  line  1,  and  x y z are  principal  axes  of  the  beam  at 
face  a,  where  xo  is  normal  to  the  beam  face.  We  define  x„ , y 
z.  as  mutually  perpendicular  unit  vectors  in  these  directions 

~3 

respectively.  Quantities  to  be  referred  to  the  local  beam  face 
coordinate  system  will  be  written  in  terms  of  these  unit  vectors. 
For  instance,  the  moment  vector  when  written  as  M = Tx  +M  v +M  z 

~ ~ 3 y ~ 3 Z ~ 3 

yields  the  torque  T about  the  local  x axis  and  the  bending 
moments  and  Mz  about  the  local  ya  and  za  axes,  respectively. 

In  Figure  2-2  (b ) the  quantities  6,  (p  , <p  are  rotations  of 

a.  y l, 

the  local  beam  face  coordinate  system  about  the  global  axes  X,  Y, 
Z.  The  angle  0 involves  beam  bending  deformation.  Later  we 
will  be  interested  in  the  increments  A0  and  A0  which  are 

y z 

incremental  bending  deformations  of  the  beam  about  the  local 
(y,z)  axes  due  to  the  moment  increments  AM  and  AM  respectively 

y z 


2 


X 


z 


Figure  2-2  a,b.  Beam  Element  Orientations 
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Associated  with  the  torque  T about  the  local  x axis  is  the  twisting 
deformation  which  is  denoted  by  ip  in  what  follows.  Also,  associated 
with  the  axial  force,  N,  acting  on  the  beam  element  is  the  elonga- 
tion 6. 

2.1  STIFFNESS  RELATIONS 

With  the  notation  just  described,  one  can  write  the  incremental 
stiffness  relation  for  the  ith  beam  element.  In  matrix  form,  this 
can  be  written  as 

A5i  = IiA?i  (2-1} 

where 

beam  element 


(2.2) 

i 

The  beam  element  stiffness  matrix,  g^,  depends  on  whether  the 
element  nodes  are  elastic  or  plastic.  The  procedure  in  Reference  1 
can  be  used  to  obtain  the  following  values  for  g^: 


AS-  are  increments  in  internal  loads  acting  on  the 
g^  is  the  element  stiffness  matrix 
AR^  are  increments  in  the  deformation  quantities 


so  that  for  element,  i,  we  have 

A0 


AR. 

~ l 


A0 

A0 

A0 

A6< 

Axp 


ay 


az 


by 

bz 


\Si 


AM 

AM, 

AM 

AM 

AN 

AT 


ay 


az 


by 

bz 
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For  elastic  behavior 


a . 


4EI 

-- 

o 

2EI 


y 


0 

4EI 

nr 

o 

2EI. 


2EI, 


4EI 


o 

2E1. 

£ 

o 

4EI. 


o 

AE 

£ 

o 


GJ 

£ 


(2.3) 


b . For  a plastic  hinge  at  the  "a"  end 


o 

o 

3EI, 

t 


o o 


o o 


3EI 


AE 

£ 


(2.4) 
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c.  For  a plastic  hinge  at  the  "b"  end 

Iei 


£ 


y 


o 

3EI 

V 


o o 


o o 


o 

AE 

£ 

o 


(2.5) 


d.  For  a plastic  hinge  at  both  ends 


li  = 


o o 


o o 


(2.6) 


AE 

£ 


We  next  define  the  G matrix  which  relates  increments  of  the 
deformation  quantities  to  the  increments  in  internal  loads  for 
all  elements  of  the  frame,  so  that 


AS  = GAR 


(2.7) 
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where,  for  a frame  having  a total  of  n^  beam  elements, 


AS 


a _ 


(2.8) 


Finally,  it  is  necessary  to  relate  the  deformation  quantities  for 
an  element  to  the  displacements  of  the  two  nodes  of  the  element, 
and  rotations  of  the  faces  at  these  two  nodes  about  the  fixed 
(global)  axes.  These  displacements  and  rotations  comprise  the 
total  degrees  of  freedom  for  each  element.  There  are  six  degrees 
of  freedom  at  each  node,  and  therefore,  twelve  degrees  of  freedom 
for  each  element.  For  n^  elements,  we  have  12  n^  equations  of 
motion.  Specifically,  for  purposes  of  solving  the  equations  in- 
crementally, we  wish  to  express  the  increments  in  the  deformation 
quantities,  AR,  to  the  increments  in  displacement  quantities,  Au, 
which  are  the  changes  in  end  point  (nodal)  locations  and  rotations. 
This  is  done  in  Section  3,  where  for  a given  element,  i,  a compat- 
ibility matrix  is  determined,  so  that  the  following  relation 
holds  for  element,  i. 


AR 


C • Au . 
~i  ~ l 


l 


Au  • 

~ l 


(2.9) 


6x12 


where 


r A 6 A0  A0,  A0,_  A 6 AiJ; 

L ay  az  by  bz 


(2.9a) 


l 
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and 


The  superscript,  T,  denotes  the  transpose  of  a matrix  quantity. 
By  appropriately  assembling  the  results  for  all  elements  of  the 
frame,  one  arrives  at  the  result 


In  Section  3,  the  C-  matrix  will  be  determined,  but  in  that  section 
the  subscript  i will  be  dropped  for  convenience,  since  all  work  in 
Section  3 refers  to  an  element. 

With  regard  to  the  equations  of  motion,  the  12  degrees  of  free 
dom  given  in  Equation  (2.9b)  will  comprise  the  dependent  variables 
for  each  element.  The  total  number  of  dependent  variables  for  the 
entire  frame  is  then  equal  to  the  dimension  of  the  vector  Au  in 
Equation  (2.10b).  We  now  wish  to  establish  the  relevant  equations 
of  motion,  which  are  to  be  solved  for  the  quantities  Au. 

2.2  EQUATIONS  OF  MOTION 

It  is  convenient  to  develop  the  equations  of  motion  by  first 
deriving  the  stiffness  matrix  for  the  associated  statics  problem 
and  then  using  the  concept  of  the  d'Alembert  force  to  account  for 
inertias.  For  the  statics  problem,  the  appropriate  stiffness 


AR  = CAu 


(2.10) 


where 


AR 


! 


(2 . 10a ,b ,c) 
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relation  to  be  used  for  our  purposes  can  be  written  as 


AP  = KAu 


(2.11) 


where  K is  defined  as  the  stiffness  matrix  of  the  entire  frame 
and  the  generalized  displacements,  Au,  are  given  by  Equation  (2.10b). 
The  vector  AP  is  the  incremental  generalized  external  loading 
vector,  so  that  at  each  node  the  concentrated  loads  correspond  to 
a particular  degree  of  freedom  in  Au.  The  stiffness  matrix,  K,  is 
now  derived  by  means  of  the  principle  of  virtual  work.  If,  for  a 
virtual  displacement  pattern  6u,  we  equate  the  internal  virtual 
work  to  the  external  virtual  work,  this  gives 


T T 

6R  S = 6u  P 


On  using  the  variational  form  of  Equation  (2.10),  this  becomes 

6uT  CT  S = 6uT  P (2.12) 

T 

For  arbitrary  Su,  this  implies  that  P = C S.  On  taking  increments, 
get 


AP  = CTAS 


(2.12a) 


Then,  on  using  Equations  (2.7)  and  (2.10)  so  that  AS  = GCAu, 
Equation  (2.12a)  becomes 

AP  = CT  GCAu  (2.13) 

When  Equation  (2.13)  is  compared  to  Equation  (2.11),  it  is  clear 
that  the  stiffness  matrix  K is  given  by 

K = CT  GC  (2.14) 


Now,  to  provide  internal  force  deformation  relations 

elasto-plastic  state  in  a typical  time  interval  x <t<x 

7 p— p+1 

internal  forces  S(t)  and  deformations  R(t)  are  written  as 
sum  of  their  value  at  x^  plus  an  additional  increment,  so 


for  an 
the 
the 
that 
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SCt)  = S(Tp)  + AS  (t ) 


(2.15) 


R(u)  = R(u(tp))  + AR(u) 


(2.16) 


where  for  an  elasto-plastic  material  the  deformation  increment,  AR, 
consists  of  two  parts  (i.e.,  elastic  Ae  and  plastic  Ar)  so 

AR  = Ae  + Ar  (2.17) 

where  the  Ar  vanish  if  the  structure  is  in  an  elastic  state. 

At  this  point,  the  incremental  equations  of  motion  are  derived 
by  starting  with  the  incremental  statics  result 

CTAS  = AP  (2.18) 

which  is  established  from  the  intermediate  virtual  work  result 

T 

given  by  Equation  (2.12a),  by  regarding  6u  arbitrary  at  that 
stage  and  going  through  the  same  argument  which  led  to  Equation 
(2.14).  In  Equation  (2.18)  we  have 

AS  = S(t)  - S(tp)  (2.19) 

AP  = P(t)  - P(x  ) (2.20) 

On  using  Equations  (2.18),  (2.19),  and  (2.20),  we  can  write 

CTS(t)  - CTS  (x  p)  = P(t)  - P(tp)  (2.21) 

In  Equation  (2.21),  since  S (x  ) and  P(t  ) ate  constants  in  the 

~ r ^ r 

interval  they  cannot  be  related  to  anything  that  is  a function  of 
time,  so  we  can  deduce  from  Equation  (2.21)  that 

P(x  ) = CT  S(x  ) (2.22) 

~ P ~ ~ P 
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Equation  (2.22)  relates  the  external  loads  at  the  pth  time  step, 

P (x  ) , to  the  compatibility  and  internal  force  matrices  C and  S 
respect  lively , at  the  t = x . On  using  Equations  (2.20)  and  (2.22), 
Equation  (2.18)  can  be  written  as 

CT  AS  - P(t)  - CT  S(x  ) (2.23) 

Now,  if  we  add  a d'Alembert  force,  -MAu  =-Mu,  (note  that  if 

u(t)  = u(x  ) + Au(t),  then  ii(t)  = Au,  since  u(x  ) is  constant)  to 

~ ~ P ~ ~ ~ P 

the  right  hand  side  of  Equation  (2.23)  the  equation  of  dynamic 

equilibrium  can  be  established.  This  equation  of  motion  can  then 

be  written  in  the  form 

MAu  + CTA  S = P(t)  - CT  S(x  ) (2.24) 

where  M is  a diagonal  lumped  mass  matrix. 

Now,  at  the  initiation  of  the  pth  time  increment  (i.e.,  at 
t = tp  ) the  matrices  G and  C are  updated  (to  account  for  elastic 
or  plastic  action  as  the  case  may  be,  and  for  geometry  changes), 
and  denoted  by  G*-P^  and  C^P^.  Then,  from  Equations  (2.7)  and  (2.10) 
we  have 

AS ( t ) = Gfp)AR(u)  - G^C^Au  (2.25) 

Next,  the  vector  Q (x  ) is  defined  by 

T 

Q(xp)  = - [C(P)]  S(Tp)  • (2.26) 

Using  updated  quantities  as  defined  in  Equations  (2.25)  and  (2.26), 
Equation  (2.24)  becomes 


MAii  + K^Au  = P(t)  + Q (x  ) 
where  the  updated  stiffness  matrix  is 


K 


tp)  _ 


G(P}  C(P} 


(2 . 27) 


(2.28) 
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3.  THREE-DIMENSIONAL  COMPATIBILITY  MATRIX 


As  discussed  in  Section  2,  the  compatibility  matrix,  C,  relates 
the  increments  in  the  local  beam  element  deformation  quantities, 

AR,  to  the  changes  in  nodal  (beam  element  end  points)  locations  and 
rotations,  Air , so  that 


AR  = CAu 


(3.1) 


where 


AR 

6x1 


C = 


:cid 


6x12 


The  elements  C—  are  determined  in  this  section.  The  notation 
indicated  in  Figure  3-1  will  be  adopted  in  what  follows. 

In  Figure  3-1,  x . y , z are  mutually  perpendicular  unit 
vectors  along  principal  axes  at  face  "a",  where  xa  is  normal  to 
face  "a"  (similar  considerations  apply  to  face  "b") . The  unit 
vector  l on  face  "a"  is  in  line  with  the  line  joining  the  nodes 
at  a and  b.  One  task  is  to  find  the  changes  in  bending  components, 
A0,ro , A0  , about  the  local  y„ , axes  due  to  small  rotations 
A<j>x,  A (py,  A <f> z of  the  local  beam  face  coordinate  system, 

(x  , y , z ) about  the  fixed  (global)  axes,  and  small  changes 

~ d ~ 3-  ~ 3. 
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in  locations  of  the  beam  element  nodes,  AX^,  AY^,  AZ^,  AX£ , AY?, 
referred  to  global  (X,Y,Z)  axes.  The  subscripts  1 and  2 refer 
to  nodes  at  face  "a”  and  "b",  respectively. 

We  will  now  show  how  the  bending  deformation  angle  increments, 

A0  and  A0  (which  are  caused  by  bending  moment  increments  AM 
ya  za  v / & ay 

and  AM  ; see  Figure  3-2)  are  related  to  A£  and  A£  which  are 
3 z a.  z a.y 

incremental  components  of  changes  in  the  £ vector. 

~ 3. 

With  regard  to  Figure  3-2,  since  x , y , z are  unit  orthogonal 
base  vectors  of  the  local  face  "a"  coordinate  system,  we  can  write 


£ = £ x + 

~a  ax  ~a 


£ y + 
ay  ca 


£ zo 
az  ~a 


(3.2) 


Then 


1 sin0  = (£  ) x (x  ) = (£  x + £ y + £ z ) x (x  ) 
a K-aJ  v~a^  v ax~a  ay  la  ~az  - a}  K~aJ 


= & r(y„  x xj  + £ (z  x x ) 
ay  a - a^  azv~a  ~a^ 


But 


y xx  = - z and  z x x = y 
la  ~ a ~ a ~ a ~ a la 


So  that  with  respect  to  the  local  system 


1 s in0 


■£  z0  + £ y 
ay  ~a  az  la 


(3.3) 


where  1 is  a unit  vector  in  the  z y plane.  Taking  the 

differential  of  both  sides  of  Equation  (3.3)  and  setting  d0, 

< 

dtay  ’ Alay’  and  d£az  = A(az  we  have 


A6a. 


cos0  A0 
a ~ a 


A£ay  ?a  + A£a  ?a 


(3.4) 


where  A0  =1  A0  This  result  can  be  written  in  matrix  form  as 

~ a ~ ~ a 

follows 
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Figure  3-1.  Face  Local  and  Global  Axes  with  Rotations 


Figure  3-2.  Beam  Element  Deformation 
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(3.5) 


Therefore,  as  was  to  be  shown,  the  deformation  angle  increments 
are  given  by  the  z and  yo  components  of  ASL  multiplied  by  the 
factor  l/cos0a.  Similar  considerations  are  used  at  node  b for 
A0^  (see  Figure  3.2),  and  in  this  case,  cos©^  would  be  involved. 

The  increments  A£  and  A£  are  defined  so  that  A£  = d£  . 

za  ya  ~a  a 

Therefore,  the  next  step  is  to  obtain  the  total  differential  of 
&a  in  terms  of  the  quantities  in  the  Au  vector.  The  reason  for 
this  is  to  enable  one  to  express  some  of  the  deformation  quantities, 
namely  A0^a  and  A0za,  terms  °f  Au,  and  in  this  manner  make  con- 

tributions to  the  compatibility  matrix  C defined  in  Equation  (3.1). 

In  order  to  obtain  d£  for  our  purposes,  it  is  convenient  to 

begin  by  defining  the  transformation  i = T£„  where  the  components 

of  are  with  respect  to  the  global  (fixed)  coordinate  system, 

while  £ is  the  unit  vector  coordinatized  in  the  a-face  frame 
~ a 

xo , y , z so  that 
a ’ 7 a ’ a 


and,  T,  which  is  the  matirx  of  direction  cosines,  is  given  by 


cos (i,x  ) 

~ ~a 

A 

cos  (3  ,x  ) 

A ~ ct 

/\ 

cos (k,x  ) 

/V  ~ d 

"Tll 

T 

12 

i 

hO 

rH 

H 

T = 

cos (i,ya) 

cos  (j  ,ya) 

cos (k,ya) 

= 

T2 1 

T 

1 2 2 

T 

23 

cos  (i’5a^1 

cos (j , za) 

cos  (k, z&) 

T3 1 

T 

1 32 

T 

33  _ 

(3.7) 
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Figure  3-3  shows  some  of  the  angles  in  the  T matrix,  where  i,  j,  k 
are  unit  vectors  in  the  global  (X,  Y,  Z)  directions.  Other 

A 

angles  not  shown  in  Figure  3-3,  such  as  (i,  y ) which  is  the  angle 
between  the  global  X axis  and  the  local  y axis,  are  defined  in  a 
similar  manner. 

Note  that  £ and  £ are  the  same  vector  described  in  different 
~ i ~ a. 

coordinate  systems.  Now,  in  terms  of  the  global  coordinates  of  the 
beam  element  nodes  (i.e.,  X^ , Y-^,  and  X2 , Y?,  Z 2 ) , one  can  ex- 
press £r  as  follows 


, 1 ? ? 7 1/9 

where  |L|  = [ (X2 -Xx ) + (Y2 - Y1 ) + (Z 2 - Z1 ) ] ' ■ 


So  that 


(3.8) 


The  constant  matrix  D is  defined  as 


-1 


o 


o 


1 


o 


o 


D 


o -1 


o 


o 


1 


o 


o 


o - 1 


o 


o 


1 


Then  £ can  be  written  as 
~ a 


~ a 


£ 


(3.9) 
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Figure  3-3.  Direction  Cosine  Angles 
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where 


X 


(3.9a) 


The  next  step  is  to  take  the  total  differential  of  £ , as  £ varies 
with  the  variations  of  the  independent  quantities  which  are  the 
elements  of  the  Au  vector.  This  means  that 

T D dX  (3.10) 

We  now  evaluate  the  three  terms  on  the  right  hand  side  of 
Equation  (3.10).  In  regard  to  the  first  term,  note  that 


represents  a scalar  change  in  length  of  the  beam  element.  Since 
L2  = (X2-Xx)2  + (Y2 -Y1) 2 + (Z  2 ‘ Z x ) 2 


d £ = 

- a 


1 

TlT 


T D X + 


1 

TlT 


dT  D X + 


Then 

2LAL=2  (X2-X1)  (AX2-AX1)  + 2 (Y2-Y1)  ( AY2  - AYX  ) + 2 ( Z 2 - Z 1 ) (AZ^AZ^ 

(3.11) 

So 


AL=  [(X2-X1)(AX2-AX1)+(Y2-Y1)(AY2-AY1)+(Z2-Z1)(AZ2-AZ1)] 

(3.12) 
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And , 


(3.12a) 


Or 


D A 


X 


(3.12b) 


where  X is  defined  in  Equation  (3.9a).  Finally,  the  first  term  in 

the  total  differential  of  i becomes 

~a 

d(p7[)T  D X = " 3 XT  DTDA  X T D X (3.13) 


— T T — • 

Since  X D DAX  is  a scalar,  Equation  (3.13)  can  be  written  as 


d(]rr)  I ? ? 


T D X XT  DT  D A X 


(3.14) 


In  order  to  evaluate  the  second  term  in  Equation  (3.10),  we  must 
determine  dT  which  indicates  how  the  transformation  matrix  T 
changes  with  rotations  of  the  face  or  local  coordinate  system  with 
respect  to  the  fixed  system.  With  regard  to  face  "a",  it  is 
shown  in  Reference  8 that  one  can  write 

AT  = T (3.15) 
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where 


0 

Aha 

-Av 

9a  = 

■S*za 

0 

4»xa 

(3.15a) 

A(J> 

_ rya 

‘Aha 

o 

The  second  term  in 

Equation 

(3.10)  is  then 

[ L|  dT  D X “ 

d T 9, 

| L | z ~ a 

D X . 

(3.16) 

This  can  be  written  in  terms  of  $ where 

~ a 

= [A<J>  A<J)  A 4>  ] 

L xa  ya  YzaJ 

So  that 

-ri-r  d T D X = -J-T-  T X*  $ 

| L | ~ ~ — | L | ~ — ~a 

Where 


X* 


o - (Z2'Z1)  ^Y2'Y1') 

(Z2-Z1)  o (X^X^ 

-(Y2-Y1)  (X2-X1)  o 


(3.17) 


(3.18) 


(3.19) 


Note  that  although  $ denotes  quantities  referred  to  face  "a", 

~ a. 

similar  considerations  hold  for  face  "b".  Finally,  the  last  term 
in  Equation  (3.10)  is  taken  as 
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(3.20) 


By  combining  the  results  in  Equations  (3.14),  (3.16),  and  (3.20) 

d£  can  be  obtained,  and  if  we  set  d£  = A£  where 
~a  ~a  ~ a 


Al 


a = 


(3.21) 


the  quantities  Al  and  Ail  can  be  found  in  terms  of  Au. 

H ay  az 

Another  component  of  the  beam  deformation  is  the  increment 
in  the  change  in  length  of  the  beam  element,  A6,  which  is  given 
by  the  quantity  AL  in  Equation  (3.11).  If  we  set  A6  = AL  and 
use  previous  notation,  we  can  write, 


1 _T  T — 

A6  = X1  D D A X 


(3.22) 


where 


The  expression  for  A6  as  given  by  Equation  (3.22)  is  in  the  desired 
form,  since  it  is  given  in  terms  of  the  changes  in  the  locations  of 
the  beam  element  nodes,  AX. 
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The  remaining  component  of  the  deformation  involves  the  twist 
angle  \p  between  the  two  faces  of  the  beam  element  as  shown  in 
Figure  3-4. 

In  order  to  compute  this  angle,  some  preliminary  considerations 
will  first  be  indicated.  If  v is  an  arbitrary  vector  then 


where  va  contains  components  of  v with  respect  to  "a"  face  local 
coordinates  (x&,  ya,  z ) , and  vr  contains  components  of  v with  re- 
spect to  the  global  X,  Y,  Z axes.  Similarly,  we  have 


v, 


Sv 


where  v^  contains  components 
coordinates  (x^,  y^ , z^)  . T 
From  Equation  (3.21)  we  have 
Equation  (3.20)  get 


of  v with 
and  S are 

rr  =~?Tyb 


respect  to  face  "b"  local 
transformation  matrices, 
so  that  on  substitution  into 


v 

~ a 


T 

T S v, 


Now,  the  angle  between  the  y-axes  (or  z-axes)  of  the  two  beam 
element  face  coordinate  systems  is  an  indication  of  the  twist 
angle  ip . Next,  select  the  y^  axis  (y-axis  on  b face)  to  be  our 
vector,  v.  Then  on  applying  Equation  (3.22)  get 


Equation  (3.23)  then  becomes 


(3.23) 


a unit  vector  by  definition. 
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Figure  3-4. 


Twist  Angle  Between  Beam  Element  End  Faces 
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Where 


(3.24) 


H 


(3.25) 


Now, 


s inf  | 


Wa  x [lb, 


(3.26) 


Where 


I f 

is  assumed  small  and  can 

X i 

yields 


la  * yb3z  la  (3.26a) 

be  neglected,  then  Equation  (3.26) 


simH  = |yba  (ya  x zfl) | 


(3.27) 


So  that  yb  is  the  magnitude  of  simp  which  is  also  equal  to 
z ^ 

or  the  (3,2)  term  of  T S (as  seen  from  Equations  3.24,  3.25  and 
3.27).  Hence , 


simp  = (3,2)  term  of  T S 


(3.28) 


However,  since  we  are  interested  in  the  increment  Aip  rather  than 

T 

simp  itself,  we  must  take  the  total  differential  of  H or  T S . 

If  we  take  the  differential  of  both  sides  of  Equation  (3.28),  we 
get 

cosipdtp  = (3,2)  term  of  d(T  S^)  . (3.29) 
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If  in  addition,  we  set  Aip  = d^  we  have 


Alp  = 3-5^7  [(3>2)  term  of  d^T  sT)].  (3. 30) 

T 

We  now  proceed  to  evaluate  dH  = d(T  S ). 
dH  = d (T  ST)  = (dT)  ST  + T d(ST) 

But , 

dT  = Tfia  and  dST  = (dS)T  = ST 

where  was  defined  in  conjunction  with  Equation  (3.11)  Therefore 

d(T  ST)  = J ST  + T flj  ST  . (3.31) 

T 

The  (3,2)  term  of  d(T  S ) is  obtained  by  performing  the  appropriate 
matrix  multiplication  indicated  on  the  right  hand  side  of  Equa- 
tion (3.31)  and  keeping  track  of  those  operations  which  contribute 

T 

to  the  (3,2)  term  of  d(T  S ).  The  result  can  be  written  in  the 
form 


(3,2)  term 
of  d(T  ST) 


where 


A6  - 


T c _t  Q 

1 32z3  1 33  z 2 
T33S2l"T31S23 
T31S22~T32S21 
T33S22~T32S23  = ~A1 


(3.32) 
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Since  we  can  now  express  all  quantities  in  the  AR  vector  in  terms 
of  those  in  Au , the  C matrix  in  Equation  (3.1),  which  is  a 6 X 12 
matrix  can  now  be  assembled  from  our  previous  results  by  adding  all 
relevant  contributions  into  their  proper  locations.  The  compatibil- 
ity matrix,  C is  denoted  as  follows: 


C11  C12  C 1 3 
r r 

c2i  ^22 


C r 

1,11  l1,12 

r 

z , 1 2 


C 


(3.33) 


r 

, 1 2 


The  elements  of  C are  given  in  the  following  section. 


3.1  ELEMENTS  OF  C MATRIX 


The  elements  of  the  beam  element  compatibility  matrix  C will 
now  be  presented.  The  following  quantities  are  first  defined 
for  convenience: 


1 

II 

X 2 

1 , - 

1 

(3.34) 

COS0  ’ 

a 

cos6^  ’ lc 

COSlfl 

X2'Xi» 

Y2 1 

Y2 " Y1 ’ Z2 1 

Z 2 ~ Z1 

(3.35) 

<X21  + 

Y2  H- 
21 

7 2 .1/2 
^21J 

(3.36) 

The  elements 


of  C are  then  given  as  follows: 
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C11  j L | 3 ( T 3 1 X 2 1 + T3  2Y21X21  + T33Z21X21 


T31Xa 


C1 2 |Ti3  (T31X21Y21  + T32Y21  + T33Z21Y21 


32  a 


C13  T7T3  (T31X21Z21  + T32 Y2 1 Z 2 1 + T33Z21 


)- 


T33Aa 


C14  " "Cll 


f = - r 
u15  l12 


C1 6 " C1 3 


'17  JEJ  ^T32Z21T33Y21') 


'18  TLT  (-T33X2l"T31Z21^) 


'19  | L I (T31Y21'T32X21^ 


Cl,10  0 


Cl,ll  0 


C1 , 1 2 0 


33 


+ 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


A 

21 

|L 

A 

22 

|L 

A 

23 

|L 

24 

" C2 1 

25 

" C2  2 

26 

~ C2  3 

A 

a 

27 

1 L | 

A 

a 

28 

|L| 

A 

a 

29 

|L| 

2,10 

= 0 

2,11 

= 0 

2,12 

= 0 

Ab 
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1 L|  3 

xb 

32 

ILI3 

T 2 1 X 2 1 + T22Y21X21  + T23Z21X21 


T21X21Y21  + T 2 2 Y 2 1 + T23Z21Y21 


T21X21Z21  + T22Y21Z21  + T23Z21 


('T23Y21'T22Z21) 

(T21Z21‘T23X21) 

^T22X21-T21Y21^ 


S31X21  + S32Y21X21 
S31X21Y21  * S 3 2 Y 2 1 


S33Z  2 1X2 1 


S33Z21Y21/ 


+ 


T21X 


T22X 


T23A 


a 


a 


a 


XbS31 

| L | 

XbS32 

I L | 


34 


Ab 

33 

1 L|  3 

34 

"C31 

35 

" C32 

36 

" C33 

37 

0 

38 

0 

39 

0 

S31X21Z21  + S32Y21Z21  + S33Z21^ 


r 

3 ,10 

. H 

|L 

c 

U3 , 1 1 

_ xb 
1 L 

C3 , 1 2 

Xb 
1 L 

X. 

C41 

1 L 

X, 

C4  2 

|L 

X- 

ii 

fO 

u 

1 L 

C44  - 

C41 

C45  ' 

" C4  2 

C46  * 

’ C43 

^S33X21"S31Z21^ 


(S31Y2i-S32X2i) 


3 ( S 2 1 X 2 1 + S22Y21X21  + S23Z21X21^  + 


3 (S21X21Y21  + S 2 2 Y 2 1 + S2 3Z 2 1X  21^ 


I 3 (S21X21Z21  + S22Y21Z21  + S23Z2l)  + ,7|3 


35 


0 


0 

0 


C 4 , 1 0 


C4 , 11 
C4 , 1 2 


^S23Y21'S22Z21^ 
yyy  (s2iz2rs23x2i^ 
jyy  ^■S22X21"S21Y21') 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


51 

52 

53 

54 

55 

56 

57 
61 

67 

68 


- x21/|l| 

- y21/|l| 

- z21/|l| 


C58  = C59  = C5,10  = C5 , 1 1 = C 5 , 1 2 
C62  = C63  = C64  = C65  C66  = 0 

XcfT32S23'T33S22) 
Xc^T33S21"T31S23^ 


0 
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C69  = Xc  (T31S22'T32S21) 

C6,10  = -C67 

C6 , 1 1 = ' C68 

C6 , 1 2 = ~ C69 


where  T — and  S — are  the  elements  of  Transformation  matrices  T, 

S defined  in  the  text. 

In  order  to  determine  A A,  and  A , we  need  cos0  , cos0K  and 

3 D C 3D 

cosijj.  They  are  calculated  as  follows: 

a.  For  cos0  , from  Equations  (3.2)  and  (3.8)  of  text, 

3 

cos0  = SL  • x = £ 
a ~a  ~a  ax 


[T11X21  + T 1 2 Y 2 1 + T 1 3 Z 2 1 


Similarly 


COS0 


[S  11^71  + S-|  7 Y 7 i + S,  -,1 


b |L|  L 11  21  12  21  1321 


b.  For  cosip,  from  Equation  (3.26a)  of  text. 


cos^  - y • (y,  ) = y = (2,2)  term  of  T S 


T 


y 
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4.  YIELD  CRITERIA  UNDER  COMBINED  LOADS 


As  the  incremental  solution  of  the  equations  of  motion 
(Eq.  2.27)  proceeds,  the  vector  Au- , for  element  i,  is  calculated 
at  each  time  step.  The  increments  in  internal  force  quantities, 

AS^,  for  element  i can  be  determined  by  using  Equations  (2.1)  and 
(2.9)  to  get 

ASptO  = g)P)c(P)  Aui(t)  (4.1) 

In  Equation  (4.1),  gT  is  the  element  stiffness  matrix  which  can 
take  on  any  of  the  forms  given  by  Equations  (2.3)  - (2.6),  depending 
on  whether  elastic  or  plastic  behavior  exists.  The  superscript,  p, 
refers  to  updated  values  at  the  pth  time  step.  The  cumulative 
values  of  the  internal  forces  of  element  i at  the  pth  time  step 
are  then  given  by 

Si(t)  = S.(tp)  + ASt(t)  (4.2) 


where 


ST  = [M  M M.  M,  N T] 
~i  1 ay  az  by  bz  J 


(4.3) 


At  a given  node  of  element  i the  bending  moments  M and  M , the 

y z 

stress  resultant,  N,  and  the  torque  T are  the  quantities  of  in- 
terest in  any  yield  cirterion,  where  they  may  take  on  critical 
values . 

For  the  loading  under  discussion,  it  should  be  possible  re- 
present the  yield  criterion  in  the  general  form 


where  Cq  is  an  appropriate  constant,  and  the  starred  quantities 
represent  fully  plastic  or  yield  values.  The  specific  form  of  the 
function,  f,  will  depend  on  the  geometry  of  the  element  cross 
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section,  the  particular  loads  which  are  assumed  to  be  chiefly 
effective  in  yielding  the  member,  and  the  kind  of  local  deforma- 
tion which  is  likely  to  occur  (this  is,  of  course,  difficult  to 
predict,  since  the  orientation  of  the  forces  cannot  be  known 
a priori).  Rather  than  attempt  to  develop  a general  criterion, 
it  appears  more  reasonable  to  derive  criteria  which  would  be 
specialized  for  the  particular  configurations  under  study,  and 
for  which  there  may  be  experimental  data  available  for  possible 
empirical  contribution  to  the  form  of  f. 

Some  simple  examples  of  possible  yield  criteria  are  presented 

here  only  for  purposes  of  giving  an  exposition  of  what  kind  of 

criterion  is  being  sought.  For  instance,  for  a rectangular 

cross  section  beam  under  combined  bending  in  one  direction 

(M  only)  and  axial  loading,  N,  bent  its  plane  of  symmetry 

(h  - beam  thickness,  b - beam  length),  the  following  yield 

g 

criterion  can  be  dervied. 


the  yield  stress.  For  the 
a method  similar  to  the 
a yield  criterion  where 
in  the  form 

(4.6) 


O ft 

where  = a bh  /4,  N = c^bh  and  a ^ is 
same  rectangular  cross  section,  by  using 
approach  in  Reference  9,  one  can  derive 

the  two  moments  M and  M are  dominant, 

l y 


= 1 


Results  for  thin  walled  sections  involving  M and  M can  also  be 

b z y 

obtained  in  a similar  manner.  Note  that  in  these  cases,  including 

* a 

that  of  Equation  (4.6),  M and  M would  need  to  be  defined 

z y 

appropriately . 

Bounding  methods  can  also  be  used  to  obtain  simple  lower 
bound  yield  loci  which  are  approximately  independent  of  cross 
section  provided  we  deal  with  thin  walled  box  sections  or  some 
types  of  I beams.  For  instance,  by  use  of  the  convexity  theorem 


39 


and  other  arguments , one  may  obtain  the  locus  in  T,  M space  which 
provides  a safe  combination  of  bending  moment,  M,  about  a "natural" 
or  "preferred"  axis,  and  torque,  T.  A pure  torsion  analysis  pro- 
vides two  points  on  the  yield  locus  T*,  o)  while  pure  bending 
analysis  gives  the  two  points  (+.  M* , o)  if  we  suppose  no  buckling 
in  compression  and  no  Bauschinger  effect.  The  fully  plastic 
moment  M*  is  proportional  to  the  yield  strength,  a . Note  that  if 
we  denote  M as  the  moment  when  the  outer  fiber  just  reaches  the 
yield  stress,  then  M*  = M^p  where  the  constant  is  a shape 
factor  which  ranges  between  the  values  of  one  and  two  in  most 
practical  situations,  and  is  closer  to  one  for  thin  walled  box 
sections.  The  four  points  (+_  T*,  o)  and  (+_  M* , o)  can  be  connected 
to  form  a quadrilateral  as  shown  in  Figure  4-1.  The  convexity 
theorem  implies  that  the  quadrilateral  locus  in  Figure  4-1  is 
"safe"  since  it  represents  a curve  closest  to  the  origin  without 
being  concave  anywhere.  The  quadrilateral,  therefore,  represents 
an  "inner"  bound  on  the  true  T,  M yield  locus,  and  applies  whatever 
the  cross  section  shape  of  the  beam. 


One  can  improve  on  this  lower  bound  in  T,  M space  by  the 
following  arguments.  For  the  torsion  problem  assume  shear 
stresses  of  magnitude  k oVer  the  cross  section  are  in  equilibrium 
with  T*;  then  lower  proportional  stresses  of  magnitude  Ak  are  in 
equilibrium  with  torque  AT*.  Similarly,  bending  stresses,  yo^  are 
in  equilibrium  with  bending  moment  yM*  (note  that  this  is  approxi- 
mate in  the  sense  that  M*  is  nearly  equal  to  M for  thin  walled 
sections).  Therefore,  if  the  combined  stresses  Ak  and  ya^  under 
say  the  von-Mises  criterion  do  not  exceed  yield,  then  the  loads 
(AT*,  yM*)  will  be  "safe".  The  von-Mises  condition  for  our  case 
gives 


(4.7) 


where  x is  the  shear  stress  and 
perpendicular  to  the  beam  axis. 


a is  a bending  stress 
Here,  x = Ak  and  a = 


on  a plane 
ya  so 

y 
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Figure  4-1.  Yield  Locus  for  Combined  Bending  and  Torsion 
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Equation  (4.7)  gives 


A2  + y2  = 1 (4.8) 

and  since  by  definition,  A = T/T*  and  y = M/M*,  Equation  (4.8) 
gives 


2 2 


The  ellipse  represented  by  Equation  (4.9)  is  shown  in  Figure  4-1, 

and  encloses  more  area  than  the  quadrilateral.  Equation  (4.9)  is 

probably  a good  approximation  for  thin  walled  box  beams,  especially 

when  the  outline  of  the  cross  section  is  closer  to  being  square. 

2 2 1/2 

The  magnitude  of  M could  be  given  by  (M  + M ) as  an  approxima- 

y z 

tion  for  general  unsymmetr ical  bending. 

In  a future  documentation,  which  is  concerned  with  computer 
assumptions  would  need  to  be  made  for  the  various  types  of  beam 
members  which  exist,  also  taking  into  account  the  fact  that  the 
response  depends  heavily  on  the  expected  loading  geometry  and 
sequence . 

In  Section  2 of  this  report,  which  is  concerned  with  computer 
implementation,  Equation  (4.9)  will  be  used  in  the  form 


For  thin  walled  sections  the  phenomenon  of  local  buckling,  or 
crippling,  is  a very  real  possibility.  How  to  take  the  gross 
effects  of  local  buckling  into  account  in  a simple  manner  is  still 
under  investigation.  As  a rough  approximation,  rather  than  use  the 
force  deformation  curve  shown  in  Figure  4-2a,  the  local  buckling 
effect  may  be  incorporated  into  the  analysis  by  use  of  the  curve 
shown  in  Figure  4-2b,  where  the  sudden  drop  in  load  is  indicated 
in  the  otherwise  simple  elasto-plastic  response.  Loading  and  un- 
loading paths  will  now  be  based  on  the  reduced  load  capacity 
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(i.e.,  lower  yield  strength).  Reduced  stiffness  caused  by  a 
smaller  effective  cross  sectional  area  may  also  be  taken  into 
account  by  smaller  loading  and  unloading  slopes,  as  shown  by 
the  dotted  lines  in  Figure  4-2b. 


FORCE 


Figure  4-2a.  Elasto - Plastic  Response  Curve 


FORCE 


Figure  4-2b.  Elasto-Plastic  Plus  Buckling  Response  Curve 
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5.  SOLUTION  OF  THE  EQUATIONS  OF  MOTION 


In  the  present  section,  a technique  for  the  numerical  time 
integration  of  the  equations  of  motion  will  be  developed.  Various 
numerical  approaches  to  the  calculation  of  the  dynamic  response  of 
complex  structural  systems  have  been  used  in  the  past  few  years. 
The  varying  degrees  of  success  of  each  method  has  been  related  to 
questions  of  accuracy  and  efficient  implementation  of  digital 
computers  (i.e.,  the  respective  algorithms  of  some  methods  lead 
to  smaller  computer  running  times  than  for  other  methods). 

Two  rather  useful  time  integration  schemes  which  have  been 

used  successfully  to  solve  complex  structural  dynamics  problems 

11  1 ? 

are  the  Houbolt  method  and  the  Newmark  Beta  method.  Houbolt's 

method  is  based  on  the  assumption  of  a cubic  curve  in  the  time 

coordinate  for  the  displacements  of  the  moving  body,  considering 

that  four  successive  ordinates  can  be  passed  through  by  a cubic 

curve.  It  is  designed  to  be  self -starting  and  unconditionally 

stable.  For  sufficiently  large  time  increments  however,  it  does 

1 3 

lead  to  some  artificial  damping  in  the  response. 

The  Newmark  Beta  method  has  been  used  successfully  in  a 
variety  of  complex  structural  response  studies,  which  involved 
treatment  of  elasto-plastic  behavior,  yield  hinges,  and  various 
dynamic  loadings  such  as  shock  or  impact,  vibration  and  earth- 
quake motion.  Some  of  the  first  applications  are  given  in 
Reference  14,  although  the  method  has  been  used  often  since 
then.  ^ The  technique  is  rather  straightforward  step  by  step 
approach,  where  the  value  of  a parameter,  3,  can  be  selected  to 
suit  the  requirements  of  the  problem  at  hand,  and  also  to  give 
unconditionally  stable  results.  The  net  effect  of  3 is  to  change 
the  form  of  the  variation  of  acceleration  in  the  time  interval. 

In  the  development  to  follow,  the  Newmark  Beta  method  will 
be  adopted.  We  shall  be  concerned  with  the  incremental  equation 
(2.27)  which  is  written  here  for  convenience.  For  the  pth  time 
step  : 
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(5.1) 


MAii  + K^P)  Au  = ptp)  + Q(-P') 


where 

q(P)  = - [CfP>]T  S(P) 


(5.2) 


Since  many  degrees  of  freedom  will  have  zero  mass  or  inertia 
associated  with  them,  Equation  (5.1)  is  partitioned  as  follows 


M A ii  + 

~a 


K 


(P) 


aa 


Au  + K 


(P) 

ag 


Aun  = P 


(P)  + 0(P)  = f(p) 
a i a 


(5.3) 


K 


(P) 

got 


Au 

-a 


+ 


= f(P) 
ig 


(5.4) 


where  the  Aua  involve  degrees  of  freedom  associated  with  masses, 
and  the  Au^  involve  degrees  of  freedom  with  no  mass  or  inertias. 
Equation  (5.4)  is  next  used  to  express  Au^  in  terms  of  Aua  as 


Au 


g 


(5.5) 


When  Equation  1 5 . 5 ) is  substituted  into  Equation  (5.3)  and  terms 
are  combined  appropriately,  the  following  equation  results. 


M ftP)  + K(P}  HP)  = f(P) 
~r  ~ 


(5.6) 


M = M 

~a 


Au 


(P) 

a 


ftP)  = 


ftp)  . 
~a 


f (P) 

~g 
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where 


-1 


K 


(P) 

33 


and  what  will  be  called  the  reduced  stiffness  matrix,  Kr , is 


K 

~r 


K 

~aa 


(5.7) 


The  appropriate  formulas  relating  the  solution  £^p^  to  £ ^P 
^(P'2)  an(j  force  quantities  can  be  derived  by  using  a technique 
described  in  Reference  15.  This  is  accomplished  here  by  first 
writing  Equation  (5.6)  for  different  time  steps  as  follows 


M + K(P)  £(p)  = f(p) 

~ ~ ~ r 


(5.8) 


M ’?(p  + 1)  + K (p+1)  ^(P  + 1)  = £(P+1)  (5.8a) 

M PP_1)  + Kr(P-l)e(p-D  = f(p-U  (5.8b) 


The  following  relations  for  velocity  and  displacement  originally 

12 

suggested  by  Newmark  are  written  next  in  terms  of  incremental 
quantities,  where  the  value  of  3 can  be  anywhere  between  zero  and 
1/4. 


pP*l)  . |( P)  + h j-pp)  + 

J(p+1)  = h j(p)  * (\  - 6) 

where  h is  the  time  increment  and 
defined  by 

£(P+1)  = U(P+1)  . U(P) 

~ ~a  ~ a 

Equations  (5.9)  and  (5.10)  can  be 


pp+l)J 

(5.9) 

h2  dP)  - Bh2  dp+1> 

(5.10) 

it  is  recalled  that  £ p+^ 


(5.10a) 

used  to  eliminate  £ ^p^  and 
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5 ^p  ^ and  obtain  the  relation, 

(P  + 1)_2  ‘^P^+  \ ('P_1  A (5.11) 

If  both  sides  of  Equation  (5.11)  are  multiplied  by  M and 
Equation  (5.8)  is  used,  one  obtains 


PP+1)  . |(p)  „ h2  ;(P)  + Bh2  U 


M 


(P+l) 


= h2f  (p) -h2K(pMP)  -8h2 
~r 


-2K(P)|(P)\ 


(5.12) 


♦ Bh2  MP+1)-K^P  + 1h(p+1)j 

■f  Bh2  /f(P"1)-K1(P'1VP"1)) 


The  terms  in  Equation  (5.12)  can  be  grouped  conveniently,  so  that 
for  the  pth  time  step,  the  following  difference  equation  can  be 
obtained 


D(p)pp)  _ b(p_ iq (p  u . Bh2dp" 2h*-p" 2Teh‘ 


(p)  + 


H 


f (p-l)+£(p-2) 


where  D^p^  = M + 8h2K  ^p^ 

~ ~ ~ X* 

B^p)  = M - (1-2  8)  h2K^ 


(5.13) 


and  the  matrix  K^p^ , has  been  calculated  at  the  end  of  the  previous 
time  step.  Note  that  the  mass  matrix,  M,  has  been  considered  con- 
stant. This  however,  is  not  necessary,  and  if  one  wishes,  one  can 
vary  M and  therefore  use  M^p).  In  fact  with  the  right  hand  side 
of  Equation  (5.13)  known,  one  regards  (5.13)  as  a system  of  equa- 
tions to  calculate  ^p-*  e Au^P^  at  the  pth  time  step.  Since  we 
are  using  a matrix  displacement  method,  the  equations  are  banded, 
and  rather  efficient  elimination  procedures  have  been  derived. ^ 
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T 

With  the  increment  [au^^J  = [Au^AuJ  known, 
vector,  u,  at  the  end  of  the  pth  time  interval  is 


the  displacement 


U(P)  = uCP'1)  + au('P')  . (5.14) 

Knowing  u^P-*,  the  results  of  Section  3 can  be  used  to  update  the 
compability  matrix  to  c(P  + ^,  for  the  next  time  step.  Also,  using 
appropriate  relations  in  Section  2 one  has 


AS(p)  = G(P)C(P)AU(P) 


where  Au^p^is  ordered  to  be  consistent  with  appropriately  defined 
matrices,  and  where  G^P^  was  calculated  in  the  previous  time  step. 
Using  the  result  in  Equation  (5.15)  the  international  forces  to  be 
used  for  the  next  time  step,  S*-P+^,  are 


S(P+1)  = S(P)  + AS(P)  (5.16) 

where  S^P-^  was  calculated  at  the  previous  time  step.  With  the 
values  of  S^,  the  yield  criteria  are  used  to  update  to  G^p  + ^^ 
for  the  next  time  step.  We  can  therefore  calculate  C^P+^  and 
g(P+1)  at-  the  encj  0f  the  pth  time  step  and  hence  can  determine 


K(P+1)  = [C(P  + 1^]  G^P  + 1^  C('P+1')  (5.17) 

In  these  conditions  the  T matrix  (Eqs.  3.7,  3.15)  must  be  updated 
as  follows: 

T(P+1)  = T(P)  + at  = + T (5.17a) 

and  similarly  for  the  S matrix.  Also,  the  philosophy  of  assembling 
K^P-*  is  discussed  at  the  end  of  this  section. 


48 


5.1  STARTING  PROCEDURE 

As  an  initial  value  problem,  the  initial  displacement,  , 

initial  velocity,  u^°-*  = and  the  time  history  of  applied 

forces  are  necessary  information.  Since  the  difference  equation, 
(5.13),  contains  terms  for  three  consecutive  time  intervals,  it  is 
required  to  express  the  displacement  u^2^  in  terms  of  initial  dis- 
placement and  initial  velocity  in  order  to  start  the  procedure. 
From  Reference  15  we  have  in  the  present  notation: 


d(»„M 


where 


E(i) 


MhC  c°)+  6h2£tl)+  /I  _g)h2ft°) 

M -U  -s)  h2^  - D(1)  - \ hV15 

p(l)_K(l)p(l)  . f(o)  = p(o) _K(o)p(o) 

~0t  ~ ~ ~f  -3 


(5.18) 


It  is  recalled  that  the  quantity  P is  a prescribed  function  and 
£^°^is  the  initial  velocity.  If  the  definition  for  and  E^^ 

are  used  together  with  the  relation  = u^’*-  u^°^  = Au^* 

° ~a  -a  ’ 

then  Equation  (5.18)  can  be  written  in  the  form 


D 


(1MD_ 


ih2>Yb(0)* 

Z ~ I ~ 


Mh£ 


r (0)  + 


3h2f 


(?  -B) 


h2f 


(o) 


(5.19) 


The  matrix  K 


(1) 

r 


is  given  by 


(5.20) 


where  the  matrices  and  represent  initial  values  to  be 

used  for  the  first  time  step.  The  quantity  e Au^0^  taken  as 

zero  since  u^  ’ itself  is  the  prescribed  initial  displacement. 
Note  also  that  = 0. 


The  solution  has  now  been  started.  The  difference  Equation 
(5.13)  can  now  be  used  to  continue  the  solution.  For  instance, 
Equation  (5.13)  can  be  used  to  obtain  writing 


49 


(5.21) 


5(2)!(2)=B(l)5a)_?(o)5(o)+  6h2  ^f(2)+  1.2  f(l)+  f(o)] 


(2)  (2)  f 2 ) 

In  order  to  calculate  E,v  , Dv  J and  are  needed  and  should  be 

calculated  before  the  start  of  the  second  time  step.  We  have 


(2)  = f(2)_K(2)£(2) 

~ a ~f  ~ B 


,(2)  = 


= M + 6 h K 


(2) 


(5.22) 

(5.22a) 


(2) 


= 


(2) 


T 


g(2)c(2) 


V. 


(5.22b) 


The  matrices  G^-2-*  and  S^-2-*  are  needed  also,  and  are  calculated 

as  follows. 

From  Equation  (5.14) 

u (1)  - A°>  ♦ Au(» 


(5.23) 


Knowing  u ^ ^ , the  compatibility  matrix  can  be  obtained  for 

the  second  time  step.  Also,  from  Equation  (5.15)  get 


AS(i)  = cCDcttW1) 


(5.24) 


Use  Equation  (5.16)  to  write 

s(2)  = s(l)  + AS(!) 


(5.25) 


where  0.  With  the  values  for  the  internal  forces  given  by 

S^,  the  yield  criteria  are  used  to  decide  on  G*'2'*  for  the  second 

12) 

time  step.  The  matrix  Kv  ' can  then  be  determined  using  Equation 

~r  r 2 ) 

(5.22c).  The  quantity,  fl  J can  also  be  determined  using 
Equations  (5.22a)  and  (5.25). 
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It  is  finally  noted  here  that  the  development  preceding  and 
leading  to  Equation  (5.17)  is  a formalism.  The  actual  details  of 
the  implementation  for  programming  purposes  will  use  the  approach 
where  the  stiffness  matrix  for  each  individual  element  is  obtained 
first,  and  then  all  elements  are  assembled  appropriately  to  obtain 
the  overall  system  matrix  The  details  will  be  given  in  a 

separate  report. 
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